
# Description -------------------------------------------------------------


# Replication Code for the Appendix to the Paper 

# "Rebels, Revenue, and Redistribution:
# The Political Geography of Post-Conflict
# Power-Sharing in Africa"

# Authors: Felix Haass, Martin Ottmann

# ---

# This script reproduces all the analytical figures and tables in the 
# appendix. For tables and figures in the main paper, see RRR_tables_plots_paper.R

# ---


# Libraries ---------------------------------------------------------------

library(lfe)
library(sf)
library(tidyverse)
library(broom)
library(stargazer)
library(countrycode)
library(haven)
library(caret)
library(doParallel)
library(geosphere)
library(Rcpp)
library(RcppArmadillo)

# Custom Functions ---------------------------------------------------------

# predict values from fixed effects linear models
source("./code/custom_functions/predict_felm.R")

# custom function to compute Conley SEs from http://darinchristensen.github.io/2015/08/30/conley/
source("./code/custom_functions/conleyses/conley.R")


# Load Data ---------------------------------------------------------------

# loading necessary data sets
load("./data/rrr_data.rda")
load("./data/rrr_data_post.rda")
load("./data/pa_only_sample.rda")
load("./data/rrr_data_full.rda")

# Data aggregated to different UoA
load("./data/rrr_adm2.rda")
load("./data/rrr_epr_growup.rda")

# Geomatched sample
load("./data/geomatched_gid.rda")


# ** TABLES ** ------------------------------------------------------------


# Table 2 -----------------------------------------------------------------


diff_in_diff_sumstats <- rrr_data %>% 
  ungroup() %>% 
  dplyr::select(nlights_1years, 
                pop_ip, 
                brd_rollsum, 
                cabinetINC,
                nonstate_rollsum) %>% 
  mutate(pop_ip = log(pop_ip + 0.01)) %>% 
  rename(`Night-Time Light Intensity` = nlights_1years, 
         `Representation in Power-Sharing` = cabinetINC,
         `Population (ln)` = pop_ip,
         `Past Battle Fatalities`= brd_rollsum,
         `Past Non-State Fatalities`= nonstate_rollsum)


stargazer(diff_in_diff_sumstats %>% as.data.frame(), digits = 2, type = "html", 
          summary = T,
          title = c("Obs", "Mean", "SD", "Min", "Max"),
          out = "./tables/appendix_table2.html", 
          float = F)


# Table 3 -----------------------------------------------------------------

rrr_data$ln_gdp_pop <- log(rrr_data$gdp_cf / 
                             rrr_data$pop_ip)



mod_gdp <- felm(log(gdp_cf / pop_ip) ~
                  log(nlights_calib_mean ) 
                | year | 0 | gid, 
                data = rrr_data %>% 
                  filter(nlights_calib_mean > 0 & gdp_cf > 0))


stargazer::stargazer(mod_gdp, type = "html",
                     out = "./tables/appendix_table3.html",
                     
                     
                     add.lines = list(c("Year FE", "\\multicolumn{1}{c}{Yes}")), 
                     covariate.labels = "log(Night Lights)", 
                     dep.var.labels = "log(GDP/pc)",
                     notes.align = "l",       
                     keep.stat = c("n", "adj.rsq"),
                     float = F,
                     align = T, 
                     dep.var.caption = NULL, 
                     notes        = "", 
                     notes.label = "", 
                     notes.append = FALSE, 
                     out.header = F)



# Table 4 -----------------------------------------------------------------

# baseline model
base_fe_time_conley <- felm(log(nlights_1years) ~
                              cabinetINC 
                            
                            | gid + year  | 0 | xcoord + ycoord, 
                            data = rrr_data , keepCX = T)

# 100km
base_conSE100 <- ConleySEs(reg = base_fe_time_conley,
                           unit = "gid",
                           time = "year",
                           lat = "xcoord", lon = "ycoord",
                           dist_fn = "Haversine", dist_cutoff = 100,
                           lag_cutoff = 5,
                           cores = 1,
                           verbose = T)
# 
# 
# 200km
base_conSE200 <- ConleySEs(reg = base_fe_time_conley,
                           unit = "gid",
                           time = "year",
                           lat = "xcoord", lon = "ycoord",
                           dist_fn = "Haversine", dist_cutoff = 200,
                           lag_cutoff = 5,
                           cores = 1,
                           verbose = T)

# 300km
base_conSE300 <- ConleySEs(reg = base_fe_time_conley,
                           unit = "gid",
                           time = "year",
                           lat = "xcoord", lon = "ycoord",
                           dist_fn = "Haversine", dist_cutoff = 300,
                           lag_cutoff = 5,
                           cores = 1,
                           verbose = T)



no_hac <- sqrt(base_conSE100$OLS)
km100 <- sqrt(base_conSE100$Spatial_HAC)
km200 <- sqrt(base_conSE200$Spatial_HAC)
km300 <- sqrt(base_conSE300$Spatial_HAC)
# km500 <- sqrt(base_conSE500$Spatial_HAC)

list_kmSEs <- list(no_hac, km100, km200, km300)


stargazer(base_fe_time_conley,
          base_fe_time_conley, 
          base_fe_time_conley, 
          base_fe_time_conley,
          se = list_kmSEs,
          type = "html",
          out = "./tables/appendix_table4.html", 
          
          covariate.labels = c("Executive Power-Sharing"),
          keep.stat = c("n", "adj.rsq"),
          add.lines = list(c("Grid-cell FE", 
                             rep("\\multicolumn{1}{c}{Yes}", 4), 
                             ""), 
                           c("Year FE", 
                             "\\multicolumn{1}{c}{Yes}", 
                             "\\multicolumn{1}{c}{Yes}", 
                             "\\multicolumn{1}{c}{Yes}",
                             "\\multicolumn{1}{c}{Yes}")),
          column.labels = c("OLS SEs", "Conley SEs / 100km", "Conley SEs / 200km", "Conley SEs / 300km"), 
          float = F,
          align = T, 
          notes.align = "l",
          dep.var.labels = "Night Lights$_{t+1}$ (log)",
          dep.var.caption = NULL, 
          notes        = "", 
          notes.label = "", 
          notes.append = FALSE, 
          out.header = F)


# Table 5 -----------------------------------------------------------------

# generate place variable
rrr_data <- rrr_data %>% 
  group_by(gid) %>% 
  mutate(excluded_treatgroup = ifelse(sum(eprINC) == 0 & 
                                        cabinetINC_treatgroup == 0, 1, 0), 
         excluded_group = ifelse(time_since_ps >= 0 & excluded_treatgroup == 1, 
                                 1, 0), 
         included_treatgroup = ifelse(sum(eprINC) >=5 & cabinetINC_treatgroup == 0, 1, 0))

# estimate baseline model
placebo_exc <- felm(log(nlights_1years) ~
                      excluded_group
                    | gid_fct + country_year | 0 | gid, 
                    data = rrr_data %>%
                      ungroup() %>%
                      filter(excluded_treatgroup == 1 | cabinetINC_treatgroup == 1) ,
                    keepX = T)

# estimate full model
placebo_exc_full <- felm(log(nlights_1years) ~
                           excluded_group +
                           log(pop_ip + 1) +
                           log(brd_rollsum + 1)  +
                           log(nonstate_rollsum +1 )
                         
                         | gid_fct + country_year | 0 | gid, 
                         data = rrr_data %>%
                           ungroup() %>%
                           filter(excluded_treatgroup == 1 | cabinetINC_treatgroup == 1) ,
                         keepX = T)


# output models to tex
stargazer(placebo_exc,
          placebo_exc_full,
          type = "html",
          out = "./tables/appendix_table5.html", 
          
          covariate.labels = c("Excluded X Power-Sharing Period", 
                               "Population (log)", 
                               "Past Battle Fatalities (log)", 
                               "Past Non-State Fatalities (log)"),
          keep.stat = c("n", "adj.rsq"),
          add.lines = list(c("Grid-cell FE", 
                             rep("\\multicolumn{1}{c}{Yes}", 2)), 
                           c("Country-Year FE", 
                             rep("\\multicolumn{1}{c}{Yes}", 2)), 
                           c("Control Group", 
                             "\\multicolumn{1}{c}{P-S Groups}", 
                             "\\multicolumn{1}{c}{P-S Groups}")),
          float = F,
          align = T, 
          notes.align = "l",
          dep.var.labels = "Night Lights$_{t+1}$ (log)",
          dep.var.caption = NULL, 
          notes        = "", 
          notes.label = "", 
          notes.append = FALSE, 
          out.header = F)


# Table 6 -----------------------------------------------------------------


paonly_model_base <- felm(log(nlights_1years) ~
                            pa_inc 
                          
                          | gid  | 0 | gid, 
                          data = pa_only_sample, keepCX = T)

paonly_model_cy_fe <- felm(log(nlights_1years) ~
                             pa_inc 
                           
                           | gid + country_year | 0 | gid, 
                           data = pa_only_sample, keepCX = T)

paonly_model_cy_fe_controls <- felm(log(nlights_1years) ~
                                      pa_inc +
                                      
                                      log(pop_ip + 1) +
                                      log(brd_rollsum + 1)  +
                                      log(nonstate_rollsum +1 )
                                    
                                    
                                    | gid + country_year | 0 | gid, 
                                    data = pa_only_sample, keepCX = T)

paonly_model_cy_fe_fullsample <- felm(log(nlights_1years) ~
                                        pa_inc +
                                        
                                        log(pop_ip + 1) +
                                        log(brd_rollsum + 1)  +
                                        log(nonstate_rollsum +1 ) +
                                        
                                        cabinetINC 
                                      
                                      | gid + country_year | 0 | gid, 
                                      data = rrr_data_full, keepCX = T)


stargazer(paonly_model_base, 
          paonly_model_cy_fe,
          paonly_model_cy_fe_controls, 
          paonly_model_cy_fe_fullsample,
          
          type = "html",
          out = "./tables/appendix_table6.html", 
          
          covariate.labels = c("Peace Agreement", 
                               "Population (log)", 
                               "Past Battle Fatalities (log)", 
                               "Past Non-State Fatalities (log)", 
                               "Representation in Power-Sharing"),
          keep.stat = c("n", "adj.rsq"),
          add.lines = list(c("Grid-cell FE", 
                             rep("\\multicolumn{1}{c}{Yes}", 4), 
                             ""), 
                           c("Country-Year FE", 
                             "", 
                             "\\multicolumn{1}{c}{Yes}", 
                             "\\multicolumn{1}{c}{Yes}", 
                             "\\multicolumn{1}{c}{Yes}",
                             "\\multicolumn{1}{c}{Yes}"), 
                           c("Sample", 
                             "\\multicolumn{1}{c}{PA Only}", 
                             "\\multicolumn{1}{c}{PA Only}", 
                             "\\multicolumn{1}{c}{PA Only}", 
                             "\\multicolumn{1}{c}{Combined}")),
          float = F,
          align = T, 
          notes.align = "l",
          dep.var.labels = "Night Lights$_{t+1}$ (log)",
          dep.var.caption = NULL, 
          notes        = "", 
          notes.label = "", 
          notes.append = FALSE, 
          out.header = F)


# ** FIGURES ** -----------------------------------------------------------


# Figure 1: Common Trends -------------------------------------------------


did_plot_line <- rrr_data %>% 
  # manually demean to remove country fixed effects
  group_by(country_year) %>%
  mutate(nlights_calib_mean = nlights_calib_mean - mean(nlights_calib_mean)) %>%
  group_by(gid) %>% 
  group_by(time_since_ps, cabinetINC_treatgroup) %>% 
  summarise(nl = mean(nlights_calib_mean, na.rm = T)) %>%
  mutate(cabinetINC_treatgroup = ifelse(cabinetINC_treatgroup == 1, "Yes", "No")) %>% 
  ungroup() %>% 
  ggplot(., aes(x = time_since_ps, 
                y = nl, 
                group = cabinetINC_treatgroup)) +
  geom_line(aes(color = cabinetINC_treatgroup, 
                linetype = cabinetINC_treatgroup), 
            size = 1.1) +
  geom_point(size = 2) +
  geom_rect(aes(xmin = -0.5, xmax = 3.5, ymin = -Inf, ymax = +Inf), 
            fill = "firebrick",
            alpha = 0.01, inherit.aes = F) +
  theme_bw() +
  scale_linetype_manual(values = c("solid", "dashed")) +
  scale_color_manual(values = c("#e31a1c", "#2b83ba" )) +
  labs(x = "Time in years before and after PS implementation",
       y = "Mean night lights\n (net of country-year effects)", 
       color='Power-Sharing', linetype = "Power-Sharing") +
  theme(legend.position = "bottom") 


ggsave(did_plot_line, filename = "./plots/appendix_figure1.pdf", 
       width = 6, height = 3.5)


# Figure 2: ADM2 + GeoEPR Maps -------------------------------------------

# code not shown


# Figure 3: ADM2 + GeoEPR Coefficient Plots -------------------------------


# * ADM2 ------------------------------------------------------------------



empty_list <- vector("list", 16)

dv_list <- c("log(sum_light + 0.01)", 
             "log(nlight_plus1 + 0.01)", 
             "log(nlight_plus2 + 0.01)", 
             "log(nlight_plus3 + 0.01)")


iv_list <- c("cabinetINC ", 
             "seniorINC", 
             "econcabINC",
             "cabinetINC_leader")

x <- 1
for(entry in iv_list) {
  for(i in 1:4) {
    fmla <- as.formula(paste(paste(dv_list[i],"~"),
                             paste(entry),
                             paste("  | ADM2ID + country_year | 0  | ADM2ID   ")))        
    model <- felm(fmla, data = rrr_adm2)
    
    empty_list[[x]] <- model
    x <- x + 1
    
    
  }
}


adm2_models <- bind_rows(lapply(empty_list, broom::tidy)) %>% 
  group_by(term) %>% 
  mutate(year = 1:n())

adm2_plot <- ggplot(adm2_models, 
                    aes(x = year, y = estimate, group = term, color = term)) +
  geom_point(position = position_dodge(width = 0.4), size = 2.5) +
  geom_errorbar(aes(ymax = estimate + 1.96 * std.error, 
                    ymin = estimate - 1.96 * std.error, 
                    linetype = term),
                size = 0.5,
                width = 0, position = position_dodge(width = 0.4)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_brewer("Type of power-sharing:", labels = c("All cabinet", 
                                                          "Leader constituency",
                                                          "Economic portfolio", 
                                                          "Senior portfolio"), 
                     palette = "Set1") +
  scale_linetype("Type of power-sharing:", labels = c("All cabinet", 
                                                      "Leader constituency",
                                                      "Economic portfolio", 
                                                      "Senior portfolio")) +
  theme_bw() +
  labs(x = "Year after begin of power-sharing", y = "Coefficient estimate") +
  theme(legend.position = "bottom")


ggsave(adm2_plot, filename = "./plots/appendix_figure3a.pdf", 
       width = 9, height = 4, scale = .9)

# * GeoEPR ----------------------------------------------------------------


empty_list <- vector("list", 16)

dv_list <- c("log(nightlight_corr_zero + 1)", 
             "log(nlight_plus1 + 1)", 
             "log(nlight_plus2 + 1)", 
             "log(nlight_plus3 + 1)")

iv_list <- c("cabinetINC", 
             "seniorINC", 
             "econcabINC",
             "cabinetINC_leader")

x <- 1
for(entry in iv_list) {
  for(i in 1:4) {
    fmla <- as.formula(paste(paste(dv_list[i],"~"),
                             paste(entry),
                             paste(" | group_fct + country_year_fct | 0 | group_fct ")))        
    model <- felm(fmla, data = rrr_epr_growup)
    
    empty_list[[x]] <- model
    x <- x + 1
    
    
  }
}


epr_sa_models <- bind_rows(lapply(empty_list, broom::tidy)) %>% 
  group_by(term) %>% 
  mutate(year = 1:n())

epr_sa_plot <- ggplot(epr_sa_models, 
                      aes(x = year, y = estimate, group = term, color = term)) +
  geom_point(position = position_dodge(width = 0.4), size = 2.5) +
  geom_errorbar(aes(ymax = estimate + 1.96 * std.error, 
                    ymin = estimate - 1.96 * std.error, 
                    linetype = term),
                size = 0.5,
                width = 0, position = position_dodge(width = 0.4)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_brewer("Type of power-sharing:", labels = c("All cabinet", 
                                                          "Leader constituency",
                                                          "Economic portfolio", 
                                                          "Senior portfolio"),
                     palette = "Set1") +
  scale_linetype("Type of power-sharing:", labels = c("All cabinet", 
                                                      "Leader constituency",
                                                      "Economic portfolio", 
                                                      "Senior portfolio")) +
  theme_bw() +
  labs(x = "Year after begin of power-sharing:", y = "Coefficient estimate") +
  theme(legend.position = "bottom")

ggsave(epr_sa_plot, filename = "./plots/appendix_figure3b.pdf", 
       width = 9, height = 4, scale = .9)



# Figure 5: Geomatching Results -------------------------------------------

diff_in_diff_geomatched <- rrr_data %>% 
  filter(gid %in% cabinc_neighbours$gid.y)


empty_list <- vector("list", 16)

dv_list <- c("log(nlights_calib_mean + 0.01)", 
             "log(nlights_1years)", 
             "log(nlights_2years)", 
             "log(nlights_3years)")

iv_list <- c("cabinetINC", 
             "seniorINC", 
             "econcabINC",
             "cabinetINC_leader")

covariates <- c("")

x <- 1
for(entry in iv_list) {
  for(i in 1:4) {
    fmla <- as.formula(paste(paste(dv_list[i],"~"),
                             paste(entry),
                             paste(" | gid + country_year | 0 | gid ")))        
    model <- felm(fmla, data = diff_in_diff_geomatched)
    
    empty_list[[x]] <- model
    x <- x + 1
    
  }
}


geomatch_models <- bind_rows(lapply(empty_list, broom::tidy)) %>% 
  group_by(term) %>% 
  mutate(year = 1:n())

geomatch_results <- ggplot(geomatch_models, 
                           aes(x = year, y = estimate, group = term, color = term)) +
  geom_point(position = position_dodge(width = 0.4), size = 2.5) +
  geom_errorbar(aes(ymax = estimate + 1.96 * std.error, 
                    ymin = estimate - 1.96 * std.error, 
                    linetype = term),
                size = 0.5,
                width = 0, position = position_dodge(width = 0.4)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_color_brewer("Type of power-sharing:", labels = c("All cabinet",
                                                          "Leader constituency",
                                                          "Economic Portfolio",
                                                          "Senior Portfolio"),
                     palette = "Set1") +
  scale_linetype("Type of power-sharing:", labels = c("All cabinet",
                                                      "Leader constituency",
                                                      "Economic Portfolio",
                                                      "Senior Portfolio")) +
  theme_bw() +
  labs(x = "Year after begin of power-sharing", y = "Coefficient estimate") +
  theme(legend.position = "bottom")

ggsave(geomatch_results, filename = "./plots/appendix_figure5.pdf", 
       width = 9, height = 4, scale = .9)


# Figure 6: Temporal Variation --------------------------------------------


# baseline + time dummies 
cabincstop_model_zero  <- felm(log(nlights_calib_mean ) ~
                                 cabINC_stop + 
                                 
                                 log(pop_ip + 1) +
                                 log(brd_rollsum + 1)  +
                                 log(nonstate_rollsum +1 )
                               
                               
                               | gid + country_year  | 0 | gid, 
                               data = rrr_data_post, keepCX = T)

cabincstop_model_1year  <- felm(log(nlights_1years ) ~
                                  cabINC_stop + 
                                  
                                  log(pop_ip + 1) +
                                  log(brd_rollsum + 1)  +
                                  log(nonstate_rollsum +1 )
                                
                                
                                | gid + country_year  | 0 | gid, 
                                data = rrr_data_post, keepCX = T)


cabincstop_model_2year  <- felm(log(nlights_2years ) ~
                                  cabINC_stop + 
                                  
                                  log(pop_ip + 1) +
                                  log(brd_rollsum + 1)  +
                                  log(nonstate_rollsum +1 )
                                
                                
                                | gid + country_year  | 0 | gid, 
                                data = rrr_data_post, keepCX = T)

cabincstop_model_3year  <- felm(log(nlights_3years ) ~
                                  cabINC_stop + 
                                  
                                  log(pop_ip + 1) +
                                  log(brd_rollsum + 1)  +
                                  log(nonstate_rollsum +1 )
                                
                                
                                | gid + country_year  | 0 | gid, 
                                data = rrr_data_post, keepCX = T)

cabincstop_model_4year  <- felm(log(nlights_4years ) ~
                                  cabINC_stop + 
                                  
                                  log(pop_ip + 1) +
                                  log(brd_rollsum + 1)  +
                                  log(nonstate_rollsum +1 )
                                
                                
                                | gid + country_year  | 0 | gid, 
                                data = rrr_data_post, keepCX = T)


cabincstop_model_5year  <- felm(log(nlights_5years ) ~
                                  cabINC_stop + 
                                  
                                  log(pop_ip + 1) +
                                  log(brd_rollsum + 1)  +
                                  log(nonstate_rollsum +1 )
                                
                                
                                | gid + country_year  | 0 | gid, 
                                data = rrr_data_post, keepCX = T)


base <- broom::tidy(cabincstop_model_zero) %>% 
  mutate(year = "+/- 0")

base_1years <- broom::tidy(cabincstop_model_1year) %>% 
  mutate(year = "+ 1")

base_2years <- broom::tidy(cabincstop_model_2year) %>% 
  mutate(year = "+ 2")

base_3years <- broom::tidy(cabincstop_model_3year) %>% 
  mutate(year = "+ 3")

base_4years <- broom::tidy(cabincstop_model_4year) %>% 
  mutate(year = "+ 4")

base_5years <- broom::tidy(cabincstop_model_5year) %>% 
  mutate(year = "+ 5")

all_time_models_cabincstop <- bind_rows(base, 
                                        base_1years, 
                                        base_2years, 
                                        base_3years, 
                                        base_4years, 
                                        base_5years) %>% 
  filter(term == "cabINC_stop")


# extract model coefficients from data
all_models <- all_time_models_cabincstop %>% 
  mutate(year = forcats::fct_relevel(year, c("+/- 0", "+ 1", "+ 2", "+ 3", "+ 4", "+ 5")))

rrr_viz_results_stop <- ggplot(all_models, 
                               aes(x = year)) + 
  geom_point(aes(y = estimate), 
             position = position_dodge(0.5)) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error, 
                    ymax = estimate + 1.96 * std.error), 
                width = 0,
                position = position_dodge(0.5)) +
  
  # scale_color_startrek() +
  geom_hline(yintercept = 0, linetype = "dashed") +
  labs(x = "Years after end of power-sharing government", y = "Coefficient Estimate") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "none")

ggsave(rrr_viz_results_stop, filename = "./plots/appendix_figure6.pdf", 
       width = 5, height = 3, scale = .9)


# Figure 7: Battle Violence Among Treated/Non-Treated Cells ---------------


brd_plot1 <- rrr_data %>% 
  ungroup() %>% 
  filter(time_since_ps < 0) %>% 
  mutate(brd_categories = cut(brd_rollsum, 
                              breaks = c(0, 1, 100, 500, 1000, 5000, 10000), 
                              right = F, include.lowest = T)) %>% 
  group_by(cabinetINC_treatgroup, brd_categories) %>% 
  tally() %>% 
  ungroup() %>% 
  complete(cabinetINC_treatgroup, brd_categories, 
           fill = list(n = 0)) %>% 
  ggplot(., 
         aes(x= as.numeric(brd_categories), y = n +1 )) + 
  geom_bar(stat = "identity", aes(fill = factor(cabinetINC_treatgroup), 
                                  group = factor(cabinetINC_treatgroup)), 
           alpha = 0.5, position = position_dodge(), 
           color = "black") +
  scale_fill_manual("", 
                    labels = c("No Rebel Constituency Areas", 
                               "Rebel Constituency Areas"),
                    
                    values = c("#252525", "#ffffff")) +
  theme_bw() +
  scale_y_log10() +
  scale_x_continuous(breaks = c(1:6),
                     labels = c("0", "1-100", "101-499", "500-999",
                                "1000-4999", "5000+")) +
  labs(x = "Battle Fatalities", y = "Number of cells (log scale + 1)") +
  theme(legend.position = "bottom") 


ggsave(brd_plot1, filename =  "./plots/appendix_figure7.pdf", 
       width = 9, height = 4, scale = 0.9)

